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1. '  INTRODUCTION 

The  direct  or  forward  problem  in  acoustics  is  the  prediction  of  the 
propagation  of  sound  based  on  specified  source,  scatterer,  and  medium 
conditions.  The  inverse  problem  is  the  deduction  of  properties  of  the 
source,  scatterer,  and  medium  from  the  propagated  field.  This  relationship 
can  be  described  with  the  mapping  of  one  function  into  another,  AX=Y, 
where  X  represents  the  specified  source,  scatterer,  and  medium  properties 
and  Y  represents  the  results  of  the  propagated  field.  The  operator  A 
provides  the  mapping  from  X  to  Y  for  the  direct  problem.  Inversion  of  A 
leads  to  the  solution  of  the  inverse  problem. 

Although  the  inverse  problem  has  been  the  subject  of  much  study, 1-4 
most  of  the  work  in  acoustics  has  been  devoted  to  the  forward  problem. 
Since  Lord  Rayleigh  predicted  the  scattered  field  from  a  sinusoidal  surface 
in  1896,  5  the  study  of  rough  surface  scattering  in  particular  has  resulted  in 
an  almost  uncountable  number  of  publications  and  about  as  many 
theoretical  developments.^”^  This  concentration  of  study  has  been  no 
less  intense  in  underwater  acoustics  where  theories  have  been  developed  to 
understand  the  scattering  from  the  ocean  surface  and  seafloor.  A  concise 
review  of  work  accomplished  in  this  area  before  1970  was  compiled  by 
Fortuin6  and  Horton.7 

The  theories  used  in  these  studies  have  been  separated  into  two  basic 
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categories  accordingto  the  mathematical  approach.  The  first  approachis 
to  model  the  scattering  with  a  distribution  of  point  scatterers,  each  with 
its  own  impulse  response  and  directivity.  Middleton8  is  regarded  as  the 
originator  of  this  rather  straightforward  approach,  often  called  the 
quasiphenomenalogical  approach.  The  second  approach,  called  the 
optical-acoustic  or  physical  approach,  was  formulated  by  Eckart.9  In  this 
method  the  scattered  field  is  described  by  Helmholtz’s  theorem  and 
evaluated  with  Kirchhoff's  boundary  conditions  resulting  in  integral 
equations  convergent  for  directional  sources. 

This  latter  approach  has  been  widely  used  by  underwater 
acousticians.  At  Applied  Research  Laboratories.  The  University  of  Texas  at 
Austin  (ARLUT),  a  group  of  investigators  have  successfully  applied  Eckart's 
theory  to  yield  theoretical  models  which  follow  experimental  data  very 
well.11"14  These  studies,  which  began  in  the  I960's,  have  produced  one 
thesis15  and  two  dissertations  16,17  at  The  University  of  Texas.  Clay  and 
Medwin 18,19  and  others20,21  have  also  verified  Eckart’s  theory  with 
experimentation. 

In  these  studies  of  the  forward  problem  the  incident  acoustic  wave 
and  the  acoustical  and  statistical  properties  of  the  scattering  medium  were 
assumed  in  an  attempt  to  predict  the  scattered  field.  The  inverse  problem 
presented  in  the  current  work  differs  in  that  a  limited  knowledge  of  both 
incident  and  scattered  pressure  is  assumed  in  efforts  to  infer  the 
scatterer’s  characteristics.  In  particular,  the  statistical  properties  of 
randomly  rough  surfaces  are  the  parameters  which  are  of  interest. 


Experimental  studies  of  the  forward  problem  at  ARUUT  were 
accomplished  with  model  surfaces  constructed  from  aeromagnetic  maps  of 
the  Canadian  Shield, 22  which  have  been  measured  for  a  particular  root  mean 
square  (rms)  surface  height  and  correlation  length.  These  surfaces  are  also 
used  in  the  inverse  scattering  problem,  the  subject  of  this  thesis.  An 
overview  of  the  acoustical  theory  supporting  the  inverse  problem  is 
discussed,  and  solutions  are  proposed  for  obtaining  estimates  of  the  rms 
height  and  correlation  length.  An  experiment  was  conducted  to  check  the 
validity  of  the  inverse  theory  and  results  from  each  of  the  surfaces  are 
compared  with  each  other  as  well  as  with  the  expected  results. 


place  on  the  data  to  reduce  the  leakage. 28 

2.  Bojarsky-Lewis  Method 

Sometimes  even  with  the  proper  choice  of  a  data  window  it  is  still 
not  practical  to  use  the  OFT  technique  to  estimate  an  FT  or  an  inverse  FT. 
For  instance,  if  the  data  record  is  so  short  that  leakage  cannot  be  corrected 
(i.e.,  through  deconvolution),  other  high  resolution  spectral  estimation 
techniques  must  be  invoked  to  estimate  the  Fourier  transform.  One  method 
is  the  Bojarsky-Lewis  method  (BL),29  which  has  been  applied  in 
three-dimensional  image  reconstruction  of  electromagnetic  wave 
scattering,  extrapolation  of  bandlimited  functions,  and  spectral 
estimation. 

The  Bojarsky-Lewis  PDF  estimate  is  derived  in  Appendix  A  and  is 
represented  with  the  summation 

J 
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.  j  =  0 . (1-1 


111.5 


where  QUO  represents  the  sampled  ve<  sion  of  Q  consisting  of  M  samples, 
*  denotes  complex  conjugate,  A£  denotes  the  even  sampling  increment  on 
Q,  and  AC=1/MA£.  The  sampled  version  of  a),  u>(jO.  will  consist  of  (1 
samples  corresponding  to  harmonically  fixed  values  of  C=jAC-  Both  Q  and 
U)  are  by  definition  periodic  functions  of  period  t1A£  and  1/A£ 
respectively.  Proper  sampling  of  Q  will  therefore  result  if  the  Nyquist 
criterion  of  1/A£>2£  max  is  followed,  where  Cmax  represents  the  maximum 

value  of  the  surface  relief  (or  a>(0=0  for  ~C.max <C^Cmax ■  a  C~>imited  (CL) 

function).  Otherwise  aliasing  of  will  occur,  and  o>(jC)  will  not 
accurately  represent  u>. 

A  major  problem  which  is  encountered  in  application  of  the  DFT  is 
the  phenomenon  of  leakage,  which  occurs  due  to  the  finite  extent  of  the 
data.  The  leakage  results  from  the  multiplication  of  Q  with  a  rectangular 
window  function  which  in  the  C  domain  corresponds  with  a  convolution  of 

a)  with  sin(2TrCmC)/TtC  (a  sinc(0  function),  where  Sm=(Smax+£min)/2- 
As  the  observation  window  Cm  approaches  infinity  the  sine  function 

approximates  a  delta  function  and  leakage  contributes  very  little  to  the 
error  in  the  DFT.  However,  it  is  not  practical  at  times  to  extend  the  data  to 
infinity,  and  so  an  optimal  window  can  be  chosen  to  place  on  the  data  to 


estimates  of  the  PDF  of  heights.  In  most  practical  situations  Q  is  not 
known  for  all  values  of  £.  A  study  for  which  the  characteristic  function  is 
known  for  various  values  of  frequency  (discrete  values  of  0  is  presented  in 
Reference  19.  The  objective  of  the  current  study  is  to  vary  the  angles  of 
incidence  and  reflection  only  over  a  modest  range  of  observation  to 
establish  measurements  of  the  characteristic  function.  So  the 
characteristic  function  is  sampled  over  a  O limited  (£0  range, 
£mjn<£<£max.  The  limits  are  selected  so  that  Q  is  bounded,  or  so  that  the 

image  solution  does  not  approach  zero  (due  to  the  finite  beamwidth  of  the 
directional  source). 

These  limitations  of  the  sampling  of  the  characteristic  function 
indicate  that  there  may  likely  be  a  resolution  problem  in  obtaining  accurate 
estimates  of  the  probability  density  function.  The  following  discussion 
outlines  the  discrete  Fourier  transform  and  two  techniques  used  in  high 
resolution  spectral  estimation  of  the  Fourier  transform. 

1.  Discrete  Fourier  Transform 

The  Fourier  transform  theory  can  then  be  extended  to  include  sampled 
functions  through  the  sampling  theorem  and  the  discrete  Fourier  transform 
(DFT).27  The  definition  of  the  DFT  pair  follows 


M-1 
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image  solution,  Eq.  11.20  becomes 


Q(0  =  <PS>/PR  Ml-* 

or 

r*  00 

0(0  =  <i>(t)  exp[-i2ffSC]  dt,  111.2 

-00 

where  the  surface  wave  number  £=2f/X  is  a  function  of  the  angles  of 
incidence  and  reflection  and  frequency.  If  the  characteristic  function  C2  is 
known  for  all  values  of  t  (all  angles  and  all  frequencies),  a),  the  probability 
density  function,  is  determined  by  the  inverse  FT 

/-+00 

0)(0=  Q(0  expli27Tttl  d£  111.3 

-  -oo 

A  sufficient  condition  for  the  validity  of  Eqs.  111.2  and  111.3  is  found 
in  Papoulis. 26  If  the  absolute  integral  of  o>(C)  is  bounded  and  finite 
(absolute  integrability),  then  its  FT  Q(0  exists  and  must  satisfy  the 
inverse  FT.  Since  the  PDF  of  surface  heights  is  bounded  for  a  finite 
surface,  its  integral  will  also  be  bounded,  and  Eqs.  111.2  and  111.3  hold. 

Thus  the  inversion  of  the  Fourier  transform  operator  will  yield 


III.  INVERSE  TECHNIQUES 


Much  can  be  Inferred  about  the  scattering  surface  from  the  two 
stochastic  quantities  derived  for  the  scattered  pressure  field.  In 
particular,  the  rms  surface  height  and  correlation  length  can  be  extracted 
from  measurements  of  the  mean  scattered  pressure  and  covariance, 
respectively.  The  rms  surface  height  is  the  square  root  of  the  mean 
squared  height  of  the  surface  profile  and  can  be  obtained  from  the 
probability  density  function  of  surface  heights.  The  correlation  length  is  a 
measure  of  the  distance  one  must  move  along  a  rough  surface  to  lose  all 
knowledge  of  one's  previous  position  and  can  be  directly  obtained  from  the 
correlation  function  of  surface  heights.  The  signal  processing  methods 
used  to  estimate  these  parameters  are  the  subject  of  this  chapter. 

A.  Probability  Density  Function 

The  probability  density  function  (PDF)  describes  the  distribution  of 
the  surface  heights  and  can  be  derived  from  the  mean  scattered  pressure 
(Section  II.B).  The  PDF  is  the  inverse  Fourier  transform  of  the 
characteristic  function;  the  conjugate  Fourier  transform  variables  are 
surface  height  t,  and  surface  wave  number  £.  In  the  surface  wave  number 
domain,  the  characteristic  function  is  defined  as  the  mean  scattered 
pressure  divided  by  the  image  solution,  PR,  of  scattering  from  a  planar 
surface.  So  the  development  is  picked  up  fromEq.  11.20  and,  fora  non-zero 


It  can  be  shown  that  for  a  suitable  random  rough  surface  (i.e., 
homogeneous,  isotropic)  the  covariance  is  dependent  only  upon  the 
difference  |r-r*|.  while  for  the  autocovariance,  the  peak  of  the  function 
occurs  at  a  zero-spatial  difference  (r=r*)  and  corresponds  with  the 
scattered  pressure  intensity.  The  covariance  decreases  rapidly  as  the  lag 
value  (spatial  difference)  increases,  the  rate  of  decrease  of  the  function 

being  greater  for  a  larger  roughness  (£rms/X  increasing).  Clay  and  Medwin 

have  shown  that  for  a  time  varying  rough  surface,  the  spatial  covariance 
reaches  a  maximum  at  a  predictable  spatial  lag  value  if  the  pressures  for 
which  the  covariance  is  calculated  are  fixed  time-delay  versions  of  the 
pressure  field.  In  other  words,  the  pressure  is  being  scattered  from  the 
same  type  of  roughness  (i.e.,  a  traveling  wave  or  surface  structure)  at 
different  times  into  different  spatial  positions.  From  calculations  of  the 
correlation  function  of  the  scattered  pressure  field,  Clay  and  Medwin  were 
able  to  predict  the  wave  velocity  of  a  traveling  surface  wave. 

For  scattering  from  a  stationary  rough  surface,  the  surface 
correlation  is  retrieved  by  spatially  correlating  the  scattered  pressure 
from  separately  insonified  areas  on  the  surface.  The  value  of  the 
correlation  function  for  a  particular  surface  displacement  is  the  peak  value 
of  the  spatial  correlation  function.  Thus  a  surface  correlation  function  can 
be  approximated  by  plotting  this  peak  value  versus  the  insonified  area 
separation. 
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assuming  et  to  be  directed  along  Rt.  The  bracketed  term  ( (...)  )  is  just  the 
pressure  field  scattered  from  a  plane  surface  corresponding  to  the  image 
solution  PR,  which  is  evaluated  with  a  stationary  phase  integral.  The  term 

<expHk#t)>  represents  the  Fourier  transform  (FT)  of  the  PDF  of  surface 
heights  or  the  characteristic  function  (CF). 

C.  Covariance  of  Scattered  Pressure 

The  covariance  of  the  pressure  field  can  be  derived  with  a 
simplification  of  Eq.  11.17.  Since  we  have  already  assumed  a  sufficiently 
directional  source,  we  can  write  Eq.  11.17  as 

Ps(r)  =  [ik/(27f)l  JJ  D0(x,y)  ({exp(ik(r*s)])/[rs])  expHktfC) 

X^xex+^yey+ez^®i dx  dy  •  1,  24 

where  r=|r|  and  s=|s|.  If  we  assume  that  the  slopes  are  negligible  the 
covariance  of  the  pressure  becomes 

<Ps(r)  Ps(r')*>  =  <exp[ik(tfC-n')l>  (k/(27rrs)]2 

XJ//J  D0(x.y)  D0(x,y)*  (ez-e,)  (e^ei')  dx  dy  dx’  dy'  11.25 

However,  the  integral  does  not  simplify  much  further  (but  can  be  evaluated 
numerically)  and  at  this  point  it  is  necessary  to  reflect  upon  the  work  of 
Reference  18  and  the  results  relating  the  covariance  to  the  correlation 
function  of  the  surface. 


quantities  we  will  be  concerned  with  are  the  so-called  mean  scattered 
pressure  and  the  covariance  of  the  pressure  field. 

B.  Mean  Scattered  Pressure 

The  ensemble  average  of  the  pressure  field,  or  mean  scattered 
pressure,  expressed  in  Eq.  11.17  is  defined  with 

<PS>  =  JJJ  Ps  ©(UxXy) «  <*h  «y  H-10 

for  to,  the  three-dimensional  probability  density  function  (PDF)  of  surface 
heights  and  slopes.  The  ensemble  average  represents  the  average  over  the 
variation  of  the  total  area  insonif led  on  the  surface.  For  a  time  varying 
surface  (C  a  function  of  time)  this  is  a  time  average,  and  for  a  fixed  surface 
(t  a  function  of  x  and  y)  this  is  an  average  over  a  number  of  insonif ied 

surface  areas,  if  C.  Cx.  and  ty  are  uncorrelated,  then  Eq.  11.17  becomes 

<Ps(r)>  =  [ik/(2TT)l  J /  D0(x.y)  ({cxp[ik(R0*Ri)l}/[R0  Ril)  <exp[-ikfftl> 
*<N;xex+Cyey+ez]-e,>  dx  dy  ,  11.19 

and  with  the  Cx  and  Cy  zero  mean  processes, 


<Ps(r)>  =  <expl-iktftl>  Uik/(2rr)l  JJ  D0(x,y)  ({exp[ik(R0*R,)]}/[R0  R,]) 


and  for  a  source  concentrating  the  incident  energy  within  a  small  region 
about  the  directional  axis  (a  directional  source  of  finite  beamwidth), 

sin  'f'o  *  sin  9j  and  sin  *  sin  ©r  11.15 

Approximations  forr0  and  rt  then  follow: 

ro  *  Ro_  C  sin  6j  and  r}  *  Rj  -  C  sin  er  11.16 

Now  the  exponential  in  Eq.  11.13  is  greatly  simplified  and  can  be 
broken  down  into  a  term  dependent  upon  the  plane  surface  geometry  and  a 
term  dependent  upon  the  surface  heights.  Also  it  should  be  noted  that  the 
denominator  r0r,  varies  slowly  with  respect  to  the  exponential  and  can  be 
replaced  with  R0Rt.  Therefore,  for  receiver  points  far  from  the  surface  and 
near  normal  incidence  (with  negligible  shadowing)  Eq.  11.13  becomes 

Ps(r)  =  [ik/(2Tr)|  / /  D0(x,y)  ({exp[ik(R0+Ri)MR0  RiD  expHkffCJ 

xlCxex<yey+ez]-e1  dx  dy  11.17 

for  #  =  sin  0j  ♦  sin  er 

Suppose  that  the  surface  heights  t(x,y)  can  be  described  with  a 
two-dimensional  zero  mean  random  variable.  This  allows  us  to  define 
various  statistical  properties  of  the  pressure  field.  The  two  particular 


(S/Sn)  (texpCikr^l/r,)  =  ik(rvei)  (CexpCikr^l/rt)  11.10 

Also,  an  expression  for  the  unit  normal  vector  and  surface  element  can  be 
derived  using  differential  geometry, 

n  =  VC/ 1  VC  |  =  {Cx e x *Cy By •* czl/ V K  )2-*<C y )2-*  11  ll.ll 

and 


d£  =  /  l(Cx)2"+(Cy)2+n  dx  dy  ,  11.12 

where  V  is  the  differential  operator,  £x  and  ty  are  the  surface  slopes  with 

respect  to  x  and  y,  and  ex,  ey,  and  e2are  the  Cartesian  coordinate  system 
unit  vectors. 

Equation  11.6  now  becomes 

Ps(r)  =  [ik/(2Tr)]//D0(x,y)  ({expIikCro+r,)]}/^  r,]) 

*KxVtyey+ezJ'*i  dx  dy  11.13 

According  to  Spetner,25  the  surface  roughness  is  assumed  smaller  than 
both  the  source  and  receiver  distances,  or 


(2C/Ro)sin  »f»0  «  I  and  (2C/Ri)sin  ^  «  1 


11.14 


another  result  is  that  multiple  scattering  and  shadowing  are  neglected. 
FromEq.  11.5,  the  Green's  function  solution  simplifies  to 

Ps(r)  =  (1/(2tt)]  JJPj(rs){(S/&n)([exp(ikr,)]/r1)}  dZ  11.6 

For  a  directional  acoustic  source  emitting  spherically  radiating 
wavefronts,  the  incident  pressure  is 

Pj  =  D0(x,y)  ([exp(ikr0)l/r0)  11.7 

in  the  farfield.  where  D0  is  the  pressure  density  proportional  to  the 
directivity  pattern  of  the  source  and  r0  the  distance  from  dl  to  the  source. 
The  partial  derivative  of  the  integrand  in  Eq.  11.6  can  be  evaluated  as 

(S/Sn)  (fexpOkr,)]/^)  =  ik(n-C|)  (Hl/(ikrt)])  ((expOkr^l/r,)  11.8 

fora  unit  normal  vector n  pointing  into  the  surface  and  et  a  unit  vector 
pointing  along  rj.  At  this  point,  the  Fresnel-Kirchhoff  approximation  can  be 
easily  applied  such  that 

kr,  » 1  ,  11.9 

which  is  equivalent  to  assuming  that  the  observation  range  is  longer  than 
the  wavelength  of  the  source.  Equation  11.8  then  becomes 


The  scattered  pressure  Ps  satisfies  the  wave  equation 

{V2  +  k2}  Ps  =  0  ,  11.2 

where  V2  is  the  Laplacian  operator  and  k=2rr/X  is  the  wave  number  for 
wavelength  X.  A  solution  to  Eq.  11.2  is 

Ps(r)  =  [I/(4tt)]  ff  Pj(rs)  {(S/Sn)  G(r,rs)}  dZ  ,  11.3 

where  Pj  is  the  incident  pressure,  (S/Sn)  the  partial  derivative  with  respect 

to  the  normal  to  the  surface  n,  dZ  the  surface  element,  and  G  the  Green's 
function  which  satisfies 

{V2  +  k2}  G  =  -4tt  S(r  -  rs)  11.4 

forG  vanishing  on  the  surface.23  The  observation  point  is  at  r  with  dZ  at 
rs.  However,  for  an  irregular  surface,  G  is  rather  hard  to  evaluate.  But  for 
a  plane  surface  placed  tangentially  to  dZ,  G  is  well  known,24  and 

{(S/Sn)  G}  =  2  {(S/Sn)  (Iexp(ikr,))/r?)}  11.5 

forr,,  the  distance  fromdZ  to  the  receiver.  This  means  that  the  radius  of 
curvature  of  the  roughness  is  assumed  much  greater  than  the  wavelength; 


II.  ACOUSTIC  THEORY 


This  section  describes  the  acoustic  theory  which  leads  to  the 
formulation  of  the  inverse  problem.  The  mathematical  framework  of  the 
forward  problem  provides  the  context  for  solution  of  the  inverse  problem. 
The  physical  model  established  by  Eckart  is  used  as  the  foundation  for  the 
inverse  theory.  The  development  used  here  follows  very  closely  that  of 
Clay  and  Medwin18  and  Boyd  and  Deaverport.14  An  expression  for  the 
pressure  field  scattered  from  rough  surfaces  is  derived  which  can  be 
siuplified  with  knowledge  of  the  incident  pressure  field  and  the 
assumption  of  a  randomly  rough  interface. 

A.  Eckart’s  Scattering  Theory 

Consider  the  scattering  geometry  of  Fig.  1;  the  source  is  spherically 
divergent,  with  position  vector  s  incident  upon  the  surface  7,  and  the 
receiver  is  at  r.  With  the  surface  height  defined  by  the  variable  t(x,y)  (in 
the  direction  of  the  z  axis)  ,  the  coordinate  system  is  oriented  such  that 
the  x-y  plane  lies  in  the  average  surface  height  described  by 

/  /  C(x,y)  dx  dy  =  0  11.1 

and  the  origin  is  the  point  where  the  beam  axis  of  the  source,  assumed 
directional,  intersects  the  mean  plane. 


v(C)  is  the  OFT  of  the  characteristic  function,  and  the  S(,)0j  are  the  prolate 

angular  wave  functions  of  the  first  kind. 

Convergence  of  the  summation  depends  upon  many  factors.  The  DFT 
is  discrete  and  thus  the  inner  product  will  be  represented  by  a  numerical 
integral  introducing  numerical  error  due  to  the  integral  approximate  The 

value  tmax  represents  an  estimate  of  the  bound  on  the  surface  heights,  and 
so  knowledge  of  the  maximum  surface  height  is  necessary.  The  SU)0j  are 

more  oscillatory  as  the  index  j  increases  and  as  the  product 
approaches  infinity.  As  a  result  most  of  the  contribution  to  the  summation 
will  be  due  to  the  lower  order  Sm0j.  Also,  if  errors  in  the  original 

characteristic  function  exist,  then  the  errors  in  the  summation  will  also  be 
large.  In  general,  Perry29  shows  that  the  BL  technique  is  numerically 
unstable  for  certain  cases.  The  degree  of  instability  depends  upon  the  value 

21^;  the  smaller  this  value  is,  then  the  more  stable  the  method  is. 

3.  Extended  Prony  Method 

Another  high  resolution  spectral  estimator  which  has  been  used  is 
the  extended  Prony  method  (EP),30  known  for  its  ability  to  accurately 
predict  the  Fourier  spectra  of  short  data  records.  The  Prony  method  has 
been  used  in  spectral  estimation,  data  reconstruction,  and  resonance 
extraction  from  transient  response  functions.  The  Prony  method,  presented 
in  detail  in  Appendix  B,  consists  of  expanding  a  complex  function  known  at 
evenly  sampled  sub-intervals  with  a  basis  set  of  complex  exponentials. 
This  expansion  can  be  expressed  as 


18 


0(0  =  1  Cj  expfejO 


111.8 


j=1 

or,  for  discrete  values  of  t  =  A  AC 

q 

Q(«)  =  l  cj[ZjlS  ,1=0 . M-l  III.9 

i=l 

which  looks  similar  to  the  DFT  expansion  of  Eq.  111.4.  However,  the  {sj }  are 
complex  in  general  and  non-harmonically  related, 

Sj  =  °<j  +  i2rrCj  111.10 


and 


Zj  =  explSjAC  ,  lll.ll 

where  the  «j  is  a  damping  factor.  Unlike  in  the  DFT,  no  assumptions  about 
periodicity  are  made  but  instead  all  parameters  are  estimated  including  the 
complex  frequencies,  S; ,  and  complex  amplitudes,  Cj.  The  frequencies  of 


the  spectra  are  not  predetermined  by  the  choice  of  data  record  length  nor  is 
the  process  restricted  to  cosines  and  sines,  if  it  is  also  assumed  that 


Cj  =  Aj  explie  j  ]  , 


are  complex  amplitudes,  then  the  FT  of  Eq.  111.25  is 


<o(0  3  2  Cj  2<xj  /l<Xj 2  +  {2tt(C~Cj  ))2)  . 


The  data  are  first  assumed  to  follow  linear  prediction  models  so 
that  a  linear  Toeplitz  system  of  equations, 


QU)  =  l  aj  Q(l-j)  ,l=q . tl-l  , 


can  be  solved  f or  the  prediction  coef  f  icients  such  that  a0=- 1.  The  equation 


£  3j  Z*<H)  =  0 


is  a  polynomial  equation  which  can  be  solved  for  the  damping  factors  and 
frequencies  Zj.  The  complex  amplitudes  are  then  retrieved  from  a  solution 

to  a  Vandermonde  system  of  equations  in  Eq.  111.9.  The  Fourier  spectrum  of 
Eq.  111.13  is  then  calculated  as  an  estimate  of  the  PDF. 

One  problem  associated  with  implementing  the  extended  Prony 
method  is  the  determination  of  the  number  of  estimation  parameters  q.  A 
method  for  doing  this  is  presented  in  Appendix  B.  Noise  also  affects  the 
accuracy  of  the  EP  method,  the  largest  impact  being  on  the  damping  factors 

{<xj}  which  become  larger  with  lower  signal-to-noise  ratios. 

B.  Correlation  Function 

The  correlation  function  can  be  obtained  from  the  covariance  of  the 
pressure  of  Section  II.C.  The  correlation  function  describes  the 
correlation  of  the  surface  profile,  i.e.,  how  much  one  portion  of  the 
surface  compares  quantitatively  with  another  portion  a  distance  away. 

However,  the  discussion  of  the  correlation  function  differs  from  the 
PDF  discussion  in  that  the  Fourier  transform  domain  consists  of  spatial 
position  and  spatial  wave  number  space,  as  opposed  to  the  surface 
heights/surface  wave  number  space  for  the  PDF. 

Following  the  definition  of  the  DFT  given  by  Eqs.  111.4  and  111.5,  the 
following  DFT  pair  involving  the  spatial  pressure  field  is  also  observed, 


M-l 

P(ftr)  =  Ap  £  P(jp)  expl-i2TTftj/M)  ,  ft  =  0 . M-l  111.16 


i 


.  j  =  0 . M-1 


111.17 


p(jp)  =  Ar 


l  P**(Hr)  expl-i2TTltj/M] 
A=0 


where  P(ftr)  represents  the  sampled  version  of  the  pressure  field  with  the 
uniform  sampling  increment  (in  Cartesian  coordinate  space)  of  Ar  (instead 
of  surf  ace  wave  number  space). 

An  estimate  of  the  covariance  of  the  pressure  field  in  Eq.  11.25  is 
possible  with  the  inverse  FT  of  the  cross-power  spectrum31  of  the 
pressure,  or 

M-1 

Rs(iU)  =  Ap  2  ss(jp)«  exp[-i2TTftj/M]  ,  It  =  0 . M-1  ,  111.18 

j=0 

where  the  covariance  is  assumed  to  be  a  function  of  the  spatial  separation 
of  the  positions  e=  |  r-r'  | .  The  cross-power  spectrum  is  defined 

ss(jp)  =  <p0(jp)  ps(jp)*>  ,  ill-19 

where  the  p's  are  the  FT's  of  the  pressure  fields  from  spots  0  and  8, 
respectively  (insonified  spots  being  separated  by  a  distance  of  8). 


An  estimate  of  the  correlation  function  of  the  surface  is  possible  by 
calculating  the  covarianceof  all  spatial  pressure  field  pairs  (p0's  andpg’s) 

associated  with  various  surface  spot  separations.  The  resulting  covariance 
can  be  normalized  with  the  autocovariance  measurements  at  £=0,  or 

M-l 

A0=Ap  l  <  |  PoOp)  |  >  111.20 

i=0 

and 

n-1 

Ag  =  ApI  <  |  P8(jp)  |  >  .  111.21 

j=0 

which  yields  the  normalized  correlation  function 

Cs(jO  =  Rs(jf)//(A0AS)  .  111.22 

The  correlation  will  have  a  peak  value  at  some  delayed  value  of  t  due  to  the 
separation  of  the  spots  5.  The  correlation  of  the  surface  can  be 
constructed  by  noting  the  peak  values  of  C  at  these  delayed  e  values  and 
plotting  them  versus  the  spot  separation  8.  This  is  verified  with  actual 
experimental  data  in  Appendix  C. 


IV.  EXPERIMENT  AND  DATA  ACQUISITION 


The  subjects  of  Chapters  II  and  111  were  the  development  of  the 
acoustical  inverse  theory  which  would  allow  inference  of  the  statistical 
properties  of  the  scattering  surface,  the  rms  height,  and  the  correlation 
length.  This  chapter  describes  the  experiment  which  is  used  to  verify  this 
theory.  First,  the  constraints  involved  in  the  theoretical  development  that 
are  most  critical  in  the  experimental  arrangement  are  outlined,  followed  by 
a  discussion  of  the  specific  experimental  parameters  used. 

A.  Considerations 

Before  the  theory  can  be  assessed,  a  valid  experiment  must  first  be 
designed  with  the  assumptions  of  the  theoretical  development  in  mind.  The 
choice  of  a  rough  surface  model  was  the  first  consideration.  Three 
pressure  release  polystyrene  models  (82  cm  by  82  cm)  constructed  from 
aeromagnetic  maps  of  the  Canadian  Shield^  were  used  and  the  statistical 
properties  of  the  surfaces  are  well  known  (see  Fig.  4).  An  actual  histogram 
of  heights  (PDF)  of  surface  3  (also  representing  surfaces  I  and  2  with 
appropriate  scale  changes)  was  developed  by  measuring  the  heights  of  1089 
points  on  the  surface  (Fig.  2).  Figure  3  represents  the  autocorrelation 
function  of  the  surface  for  two  perpendicular  orientations  of  the  surface. 
A  Gaussian  and  exponential  curve  have  been  fit  to  both  the  histogram  and 
correlation  function,  respectively.  A  measure  of  the  rms  height  is  obtained 
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SURFACE  HEIGHTS  (cm) 


I 
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FIGURE  2 

HISTOGRAM  OF  HEIGHTS  MEASURED  AT  1089  POINTS  IN  THE  CENTRAL 
QUARTER  OF  THE  THIRD  ROUGH  SURFACE  (COURTESY  OF  S.  K.  MITCHELL) 

STANDARD  DEVIATION  OF  THE  GAUSSIAN  CURVE  IS  1.13  cm,  WHICH 
APPROXIMATES  THE  rms  HEIGHT  OF  THE  SURFACE 
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FIGURE  3 

AUTOCORRELATION  OF  THE  SURFACE  MEASURED  FROM  THE  CENTRAL 
QUARTER  OF  THE  THIRD  MODEL  SURFACE  (COURTESY  OF  S.  K.  MITCHELL) 

THE  E  FOLDING  VALUE  OF  THE  EXPONENTIAL  CURVE  IS  4.51  cm,  WHICH 
APPROXIMATES  THE  CORRELATION  LENGTH  OF  THE  SURFACE 


AS-85-309 


with  the  standard  deviation  of  the  Gaussian  function,  and  a  measure  of  the 
correlation  length  is  obtained  by  observing  the  e-folding  (1/e)  point  of  the 
exponential  function. 

The  next  consideration  was  that  of  establishing  an  acoustic 
measurement  geometry.  Recall  that  the  acoustic  source  is  restricted  to 
near-normal  incidence  so  that  there  will  be  a  minimum  shadowing  of  the 
surface.  Shadowing  studies32  were  conducted  on  the  first  three  rough 
surfaces  and  estimates  of  the  shadowing  functions  obtained.  A  choice  of 

9j>50°,  such  that  the  shadowing  function  is  unity,  will  assure  that 

shadowing  is  insignificant.  Figure  4  represents  the  experimental  geometry 
and  is  referenced  for  the  rest  of  the  discussion. 

The  wavelength  of  the  incident  acoustic  pressure  must  also  be 
selected.  A  rough  guide  to  this  choice  is  the  Rayleigh  criterion10  of 
surf  ace  roughness.  A  surf  ace  is  considered  ef  f  ectively  smooth  if 

Ws/X  0  or  9j  -»  0  ,  IV.1 

which  should  set  an  upper  limit  on  X.  Another  guide  to  the  choice  is  to 
select  a  range  on  £  for  sampling  the  main  structure  of  the  characteristic 
function.  If  a  Gaussian  PDF  is  assumed,  the  use  of  the  3Crms  point  of  the 
characteristic  function  yields  a  lower  limit  onX,  or 

3/Crms  >  (sinei  +  sinM/x  IV.2 
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(xp,yp,zp)  (xr.yr!  to  yr2,zr) 


Model 

1 

2 

3 

rms  Height,  Crms 

0.231  cm 

0.462  cm 

0.924  cm 

Correlation  length 

6.48  cm 

6.48  cm 

6.48  cm 

Wavelength,  X 

1.88  cm 

1.88  cm 

1.88  cm 

Transmit,  (xp,yp,zp) 

(0,-50,275) 

(0,-50,275) 

(2,-70,293) 

Receive,  (xr,yrl;yr2,zr) 

(0,27;73,275) 

(0,27,73,275) 

(-3t-l.12l.273) 

Spot  radius.  a«b 

16.5  cm 

16.5  cm 

16.5  cm 

Number  of  spots 

49 

49 

69 

Transmit  pulsewidth 

450  ps 

450  ps 

600  ps 

Repetition  rate 

25  /s 

25  /s 

30  /s 

Receive  pulsewidth 

125  ps 

125  ps 

200  ps 

FIGURE  4 

EXPERIMENTAL  GEOMETRY  AND  PARAMETERS 
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The  solutions  are  valid  for  directional  sources  only  so  a  fairly  small 
beamwidth  is  necessary.  This  parameter  is  of  particular  interest  since  the 
beamwidth  along  with  the  range  will  determine  the  illumination  spot  size 
on  the  surface  and  the  number  of  independent  spots  which  can  be  insonified. 
The  larger  the  number  of  spots,  the  larger  will  be  the  ensemble  over  which 
the  mean  scattered  pressure  is  calculated  in  Eq.  11.20.  Also,  the  dimension 
of  the  spots  should  be  on  the  order  of  a  correlation  length  so  that  at  least 
one  correlation  length  of  the  surface  is  insonified. 

The  insonified  area  on  the  surface  is  considered  to  be  determined  by 
the  -3  dB  beamwidth  (0)  contour.  For  a  directional  source  beam  pattern 
which  exhibits  azimuthal  symmetry,  the  spot  is  elliptical  with  semi-major 
axis  a  and  semi-minor  axis  b, 

a=  zp  tan(0/2)/sin2e  j  IV.3 

and 

b=asin0j  IV.4 

■  i 

where  zp  is  the  range,  being  determined  by  the  farfield  distance  of  the 
source,  or  zp>S/X  for  S,  the  acoustically  active  surface  area  of  the  source. 

The  orientation  of  the  source  and  receiver  should  also  be  such  that 
the  source  does  not  interfere  with  the  measurements  of  the  field.  So 
restricting  the  measurements  to  the  -10  dB  contour  will  yield  an  elliptical 


spot  of  2/ (10/3)  times  that  of  the  surface  spot.  This  will  also  prevent 
the  characteristic  function  from  being  undefined  as  the  image  solution 
approaches  zero. 

Also,  a  time  Tj  is  required  to  illuminate  the  surface  within  the  -3  dB 
contour.  From  Fig.  4  it  is  noted  that 

R„  =  yp/sin(©  j+^/2) 

R2,  =  yp/sinOj-£/2) 

and 

Tj  =  (R2rRii)/c  ,  IV.7 

where  c  is  the  sound  speed  in  water  (1500  m/sec).  This  time  delay  results 
in  a  rise  time  that  will  allow  a  pulsed  signal  to  be  detected  as  a  steady 
state  signal.  There  is  also  a  time  delay  Td  between  transmit  at  the  source 
and  receive  at  the  receiver,  or 

T(j  =  (R2j4Rh)/c  IV. 8 

The  experimental  measurements  were  done  in  a  sonar  model  tank 
room,  which  meant  that  there  would  be  multiple  reflections  due  to  the 
walls,  water  surface,  and  tank  bottom.  This  problem  was  overcome  by 


IV.5 

IV.6 


operating  the  source  in  a  pulsed  mode  rather  than  a  continuous  mode  so 
that  a  particular  pulse  rate  and  pulse  width  can  be  selected  such  that  the 
reflections  will  not  interfere  with  the  return  pulse  from  the  rough  surface. 
Thus,  the  timing  involved  in  obtaining  the  scattered  pressure  measurements 
is  critical  in  that  the  proper  delay  of  the  observation  interval  should  be 
chosen  to  capture  the  valid  scattering. 

It  is  impractical  (if  not  impossible)  to  experimentally  measure  the 
characteristic  function  for  all  possible  surface  wave  numbers.  However, 
we  can  obtain  a  finite  number  of  discrete  values.  According  to  the 
sampling  theorem  he  characteristic  function  must  be  sampled  at  least  at 
the  Nyquist  rate  in  even  intervals  of  the  variable  £.  Thus  the  receiver  and 
source  positions  or  wavelengths  must  be  selected  such  that  C  is 
incremented  in  even  intervals.  This  is  done  by  maintaining  the  source  as 

stationary  (6j  fixed),  radiating  one  frequency  (X  fixed),  and  moving  the 

receiver  (er  varied)  such  that  Q  is  sampled  in  even  increments.  The 

covariance  function  requires  the  sampling  of  the  pressure  field  in  even 
increments  of  space,  so  a  technique  must  be  developed  to  allow  sampling 
both  in  space  and  surface  wave  number  space,  £. 

B.  Geometry  and  Experimental  Equipment 

The  forward  scattered  data  were  collected  in  the  sonar  model  tank 
room  at  ARL:UT.  The  scattering  surfaces  were  a  plane  pressure  release 
surface  and  three  randomly  rough  surfaces  constructed  specifically  for 
scattering  studies.  The  complex  pressure  field  was  measured  using  a  line 
and  cone  transducer33  as  the  narrowbeam  projector  and  an  H-56  standard 
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rms  Height  (in  cm) 


Surface 

Actual 

DFT 

Elementary 

Bandwidth 

BLtl 

EPM 

I 

.231 

219 

214 

71.4 

32.7 

2 

.462 

237 

214 

82.5 

41.2 

3 

.924 

47.0 

22.4 

8.81 

1.40 

Correlation  Length  (in  cm) 

Surface 

Actual 

Calculated  Spatial 

Spot  Separation 

Lag 

Increment 

1 

6.48 

7.03  .55 

3.18 

2 

6.48 

3.32  .55 

3.18 

3 

6.48 

5.18  .55 

2.29 

TABLE!:  Summary  of  results  of  the  estimates  of  rms  height  and 
correlation  length 


the  Gaussian  curve,  carefully  noting  the  degree  of  goodness  of  fit  between 
the  estimates  and  analytical  curve.  For  more  realistic  comparisons, 
however,  the  structure  about  the  central  maximum  must  also  be  fitted  with 
a  Gaussian  curve. 

The  accuracy  of  the  rms  height  estimates  is  a  problem  due  to  the 
finite  extent  of  the  observation  window.  The  elementary  bandwidth  is  the 
inverse  of  the  observation  window  and  is  directly  related  to  the  resolution 
of  the  Fourier  transform  estimate.  Thus,  the  length  of  the  data 
observation  window  will  affect  each  of  the  three  techniques  and  resolution 
will  improve  as  the  observation  window  is  increased.  Also,  greater 
resolution  is  achieved  with  EPM  than  with  BLN,  and  better  resolution  with 
BLN  than  DFT.  All  three  methods  are  somewhat  sensitive  to  changes  in  the 
rms  height  when  the  observation  window  is  the  same,  as  seen  in 
comparisons  of  the  PDF  estimates  of  surfaces  1  and  2.  These  results  are 
summarized  in  Table  I,  the  values  representing  the  estimated  rms  heights 
obtained  from  Gaussian  curves  fitted  to  the  PDF  estimates  of  the  three 
techniques. 

1.  DFT  Technique 

Figures  9,  10,  and  11  represent  the  discrete  Fourier  transform 
technique  applied  to  the  characteristic  functions  of  surfaces  1,  2,  and  3, 
respectively.  These  results  will  be  the  benchmark  by  which  the  other  two 
techniques  are  compared  for  resolution.  For  surfaces  1  and  2,  the 
elementary  bandwidth  of  the  surface  heights  is  214  cm  when  the  observed 
data  record  is  used  by  itself.  However,  this  resolution  of  heights  can  be 


42 


least  means  square  solution  for  the  {Cj }  in  Eq.  III.9.  The  Fourier  spectrum 

is  then  computed  with  application  of  Eq.  111.13,  and  as  in  the 

Bojarsky-Lewis  method,  extrapolation  or  interpolation  of  the 

characteristic  function  is  possible. 

The  correlation  function  computation  proceeds  with  the  FFTs  of  the 
scattered  pressure  fields  from  each  of  the  spots  (Eq.  111.17).  Since  the 
correlation  function  estimate  is  valid  for  zero-mean  stochastic 
processes,31  the  mean  value  of  the  pressure  must  be  zero  before 
proceeding.  The  auto-power  and  cross-power  spectra  of  the  pressure 
fields  corresponding  with  spots  separated  a  given  uistance  are  then 
computed  with  Eq.  Hl.19,  the  ensemble  average  being  the  average  over  the 
ensemble  of  spot  pairs.  The  normalized  correlation  function  (Eq.  111.22)  for 
each  possible  spot  separation  is  computed  via  the  inverse  FFT  and  with 
knowledge  of  the  zero-lag  values  of  the  auto-power  spectra  (Eqs.  IH.20  and 
111.21).  The  peak  value  of  the  correlation  is  then  plotted  versus  the  spot 
separation  of  the  spots. 

B.  PDF  Estimates 

This  section  presents  the  numerical  solutions  resulting  from  the 
application  of  the  three  PDF  estimation  techniques  upon  the  experimental 
data.  The  results  are  presented  for  each  of  the  three  model  rough  surfaces. 
The  general  structure  of  the  PDFs  is  bounded  and  Gaussian  in  shape. 
Therefore,  a  Gaussian  function  is  fit  to  the  PDF  estimates  to  minimize  the 
mean  squared  error  between  the  Gaussian  and  the  estimate  (a  least  mean 
square  fit).  The  rms  height  value  is  taken  to  be  the  standard  deviation  of 
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the  fast  Fourier  transform  (FFT)  algorithm. 27 

For  the  application  of  the  Bojarsky-Lewis  method,  the  prolate 
angular  wave  functions  must  be  calculated  for  use  in  Eq.  111.6.  These 
functions  are  calculated  with  knowledge  of  and  Cmax  and  the  algorithm 

of  Reference  36.  The  value  used  for  is  the  average  value  of  the  surface 
wave  number  over  the  range  of  sampling  of  the  characteristic  function,  and 
the  value  used  for  Cmax  is  an  estimate  of  the  bound  on  the  surface  heights. 

This  value  is  assumed  to  be  three  times  the  actual  rms  surface  height  of 
the  model  surface.  Recalling  the  truncated  summation  of  Eq.  111.6,  one 
notes  the  inner  product  between  the  prolate  angular  wave  functions  and  the 
estimate  of  the  PDF.  v(0  of  Eq.  111.7.  The  estimate  of  v  is  computed  with 
the  FFT  of  Eq.  111.5  as  in  the  discrete  Fourier  transform  technique  above 
(however,  no  optimal  window  is  applied).  The  inner  product  is  then 
calculated  as  a  numerical  integral.  Due  to  the  continuous  nature  of  the 
estimate  of  Eq.  111.6,  any  interpolated  or  extrapolated  value  for  the  PDF  is 
possible,  the  only  limitation  being  the  truncation  of  the  summation. 

As  in  the  Bojarsky-Lewis  method,  the  extended  Prony  method  is  a 
multiple  step  process.  The  sampled  characteristic  function  forms  the  data 
vectors  of  Eq.  B.ll  for  calculation  of  the  data  matrix  of  Eq.  B.17,  thus 
requiring  knowledge  of  the  value  forq.  The  eigenvalue  decomposition  of 
the  matrix  is  computed  and  the  weakest  eigenvector  (corresponding  to  the 
smallest  eigenvalue)  is  the  vector  of  Eq.  B.I8.  The  polynomial  of  Eq.  111.15 

is  then  rooted  to  find  the  values  for  the  (Zj),  which  are  then  used  to  find  a 


pressure  corresponds  to  the  -3  dB  beamwidth  of  the  image  solution.  The 
scattering  from  surface  3  is  observed  over  the  -10  dB  beamwidth  of  the 
image  solutioa  This  data  observation  length  determines  the  resolution  of 
the  PDF  estimates,  as  will  be  seen. 

The  plots  were  digitized  so  that  data  points  could  be  interpolated 
for  sampling  the  pressure  in  even  increments  of  both  surface  wave  number 
space  (0  and  Cartesian  coordinate  space  (r).  A  “cubic  spline"  technique 
(fitting  a  cubic  equation  to  data  points)  was  used  for  interpolation  of  the 
digitized  amplitude  and  a  linear  interpolation  method  was  used  for  the 
digitized  phase. 

The  uniform  sampling  in  surface  wave  number  space  allows  the 
calculation  of  the  PDF  and  the  sampling  in  coordinate  space  allows  the 
calculation  of  the  correlation  function.  Thus,  following  interpolation, 
solutions  to  the  inverse  problem  can  be  sought.  For  the  calculation  of  the 
PDF,  the  characteristic  function  must  be  computed  using  Eq.  111.1.  The 
numerator  is  the  mean  scattered  pressure  and  is  calculated  with  the 
ensemble  average  of  the  pressure  field  over  the  independently  insonified 
spots.  The  image  solution  is  the  pressure  field  for  the  reflection  from  the 
planar  surface.  The  inverse  Fourier  transform  of  the  characteristic 
function  must  now  be  estimated  with  one  of  the  three  techniques  of 
Chapter  111. 

Before  the  application  of  the  discrete  Fourier  transform,  an  optimal 
window,  the  Kaiser-Bessel  window, \s  applied  to  the  characteristic 
functioa  The  discrete  Fourier  transform  of  Eq.  111.5  is  then  applied  using 


FIGURE  8 

FLOWCHART  OF  THE  PROCESSING  OF  THE  SCATTERED  PRESSURE 
DATA  FOR  INFERENCE  OF  THE  PDF  AND  CORRELATION  FUNCTION 


V.  RESULTS 


This  chapter  describes  the  analysis  of  the  scattered  pressure  field 
measurements  using  the  processing  techniques  discussed  in  Chapter  III.  The 
numerical  results  of  the  calculations  for  each  of  the  model  rough  surfaces 
are  presented  and  compared  with  the  actual  statistical  parameters  of  the 
model  surfaces.  Suggestions  are  made  for  improving  the  accuracy  of  the 
results  where  appropriate. 

A.  Data  Processing 

The  implementation  of  the  signal  processing  techniques  in  Chapter  HI 
is  straightforward,  in  that  the  equations  necessary  for  the  processing  are 
all  present  in  the  text.  Figure  8  is  a  flowchart  of  these  processing 
techniques,  which  were  programmed  in  FORTRAN  for  use  on  a  CYBER  171 
digital  computer  available  at  ARLDT.  Although  the  actual  coding  is  not 
presented  here,  archived  copies  of  the  programs  along  with  documentation 
for  usage  are  available  upon  request. 

The  complex  pressure  field  measurements  are  represented  with  the 
amplitude  and  phase  sweeps  obtained  from  the  receiving  hydrophone  (the 
output  of  the  voltmeter  and  phasemeter35).  The  scattered  pressure  field 
from  both  the  planar  surface  and  the  model  surface  are  necessary  to 
compute  the  PDF,  the  scattering  from  the  plane  surface  being  used  as  the 
image  solution.  For  surfaces  1  and  2,  the  observation  interval  of  the 


FIGURE  7 

A  SAMPLE  PHASE  PLOT  OF  THE  OUTPUT  OF  THE  PHASEMETER  versus  THE 
SPATIAL  POSITION  OF  THE  HYDROPHONE 

THE  SAMPLING  INCREMENT  IS  UNIFORM  IN  CARTESIAN  COORDINATE  SPACE 
THIS  REPRESENTS  THE  PHASE  OF  THE  PRESSURE  FIELD  OF  THE  SCATTERING 
FROM  ONE  INSONIFIED  SURFACE  SPOT  ON  SURFACE  THREE 
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AMPLITUDE 


SPATIAL  POSITION  (cm) 


FIGURE  6 

A  SAMPLE  AMPLITUDE  PLOT  OF  THE  OUTPUT  OF  THE  VOLTMETER  versus  THE 
SPATIAL  POSITION  OF  THE  HYDROPHONE 

THE  SAMPLING  INCREMENT  IS  UNIFORM  IN  CARTESIAN  COORDINATE  SPACE 
THIS  REPRESENTS  THE  AMPLITUDE  OF  THE  PRESSURE  FIELD  OF  THE  SCATTERING 
FROM  ONE  INSONIFIED  SURFACE  SPOT  ON  SURFACE  THREE 


LINE-AND-CONE  STANDARD  H-56 


FIGURE  5 

DATA  ACQUISITION  EQUIPMENT  FOR  EXPERIMENTAL  MEASUREMENTS 
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from  the  receiver  were  amplified  and  filtered  for  analysis  of  phase  and 
amplitude  using  a  pulse  phasemeter34  and  voltmeter.  The  amplitude  and 
phase  were  recorded  using  chart  recorders  in  the  sweep  mode,  with  the 
sweep  synchronized  with  the  scanning  of  the  field.  Figure  5  is  a  block 
diagram  of  the  equipment  used  for  data  acquisition.  The  continuous 
amplitude  and  phase  plots  resulting  from  the  above  experiments  were  then 
digitized  for  data  storage  and  subsequent  processing,  examples  of  which 
are  shown  in  Figs.  6  and  7.  Figure  6  is  an  amplitude  plot  of  the  scattered 
pressure  from  one  surface  spot  insonified  on  surface  3.  Figure  7  is  the 
phase  plot  corresponding  to  the  amplitude  plot  of  Fig.  6. 


33 


hydrophone  as  the  omnidirectional  receiver.  The  projector  remained  fixed 
at  a  given  angle  of  incidence  upon  the  surface  and  the  receiver  was  moved 
to  form  an  array  of  pressure  measurements.  The  pressure  field  was 
scanned  for  each  insonified  area  on  the  model  surface  to  form  an  ensemble 
of  pressure  fields.  The  projector  and  receiver  coordinates  were  selected 
so  that  shadowing  was  minimized.  When  a  spot  size  on  the  order  of  a 
correlation  length  was  insonified  neither  the  projector  nor  receiver 
interfered  with  the  other  (see  Fig.  4). 

Due  to  the  finite  size  of  the  water  tank,  a  pulsed  cw  signal  was 
transmitted  at  a  frequency  of  80  kHz.  The  pulse  was  rectangular  with  a 
fixed  width  and  repetition  rate.  The  projector  had  a  -3  dB  beamwidth 
covering  a  spot  of  radius  16.5  cm  on  the  surface  and  an  acoustic 
wavelength  of  1.88  cm. 

The  model  surface  was  located  in  the  farfield  of  the  projector  at 
near-normal  incidence  producing  a  -3  dB  insonification  spot  size 
approximating  a  circle.  A  number  of  spots  were  insonified  leaving  an  area 
the  size  of  one  spot  radius  uninsonified  at  the  edge  of  the  surface.  The 
total  area  insonified  was  more  than  half  the  surface  so  that  the  results 
from  the  processing  could  be  compared  with  the  physical  characteristics  of 
Figs.  2  and  3. 

The  receiver's  active  acoustic  size  was  approximately  one 
wavelength  in  dimension.  The  receiver  was  moved  by  a  motorized  column 
in  a  direction  parallel  with  the  model  surface.  Only  the  steady  state 
portion  of  the  scattered  pulse  was  sampled  for  processing.  The  outputs 
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improved  by  assuming  the  characteristic  function  is  zero  outside  the 
observation  window  and  computing  the  FFT  with  the  extended  data  record 
thus  improving  the  elementary  bandwidth.  Although  this  is  cheating,  the 
resolution  of  the  estimate  is  not  improved  because  of  the  leakage  due  to 
the  finite  length  of  the  non- zero  data.  The  only  improvement  is  in  the 
view  of  the  finer  structure  of  the  PDF. 

However,  an  improvement  in  the  resolution  can  be  seen  in  the  PDF 
estimate  of  surface  3  for  which  the  elementary  bandwidth  is  22  cm.  The 
closest  estimate  of  the  rms  height  (from  the  DFT)  is  represented  with  the 
Gaussian  curve  fit  of  Fig.  9  and  is  47.0  cm  compared  to  the  true  rms  height 
of  0.924  cm.  It  is  obvious  that  high  resolution  techniques  are  necessary 
for  better  estimates  of  the  rms  height. 

2.  BLfl  Technique 

Figures  12,  13,  and  14  reflect  the  application  of  the  Bojarsky-Lewis 
method  to  the  characteristic  functions  of  surfaces  1,  2,  and  3.  It  should  be 
noted  that  the  PDF  estimates  of  Figs.  12  and  13  no  longer  fit  a  Gaussian 
shape  due  to  the  appearance  of  secondary  structure.  The  rms  estimates 
with  the  least  mean  square  Gaussian  fit  are  therefore  deceptive,  since  it  is 
the  central  structure  of  the  PDF  which  should  be  observed.  For  a  Gaussian 
curve  fit  to  the  central  portion,  rms  height  estimates  of  71.4  and  82.5  cm 
are  obtained  for  surfaces  1  and  2.  Although  the  Bojarsky-Lewis  method 
makes  use  of  the  discrete  Fourier  transform,  an  improvement  in  resolution 
is  evident,  especially  in  the  PDF  for  surface  3.  The  rms  height  for  this 
surface  is  8.81  cm,  a  definite  improvement  in  the  estimate.  However,  it 
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FIGURE  12 

THE  PROBABILITY  DENSITY  FUNCTION  FOR  THE  FIRST  ROUGH 
SURFACE  VIA  THE  BOJARSKY-LEWIS  METHOD 

STANDARD  DEVIATION  OF  THE  GAUSSIAN  CURVE  IS  224  cm 
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FIGURE  13 

THE  PROBABILITY  DENSITY  FUNCTION  FOR  THE  SECOND  ROUGH 
SURFACE  VIA  THE  BOJARSKY-LEWIS  METHOD 

STANDARD  DEVIATION  OF  THE  GAUSSIAN  CURVE  IS  208  cm 
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still  does  not  compare  with  the  true  value  of  0.924  cm.  So  it  is  necessary 
to  use  a  higher  resolution  technique. 

3.  EPM  Technique 

Figures  15,  16,  and  17  show  the  extended  Prong  technique  applied  to 
the  characteristic  functions  of  surfaces  1,  2,  and  3.  When  applied  to 
surfaces  1  and  2,  it  is  noted  that  the  least  mean  square  Gaussian  curve  does 
not  fit  well.  When  a  Gaussian  curve  is  fit  to  the  central  portion  of  the  PDF 
visually,  estimates  of  32.7  and  41.2  cm  are  obtained  for  surfaces  1  and  2. 
As  is  expected  with  the  longest  data  window,  the  best  rms  height  estimate 
of  3.2  cm  is  obtained  for  the  third  surface.  Fitting  the  best  Gaussian  curve 
to  the  central  portion  yields  an  estimate  of  1.4  cm.  The  steps  necessary 
for  further  improvement  of  the  estimates  are  obvious  when  comparing 
surfaces  1  and  2  results  to  surface  3  results.  The  extension  of  the  length 
of  the  data  record  improves  the  resolution  as  expected.  Therefore,  if  the 
sampling  of  the  characteristic  function  is  extended  past  the  -10  dB  points 
on  the  image  solution,  a  further  improvement  in  the  estimates  of  the  PDF 
should  result.  The  only  other  alternative  is  to  use  a  higher  resolution 
technique  for  estimation  of  the  PDF. 

C.  Correlation  Estimates 

The  surface  correlation  functions  were  constructed  from  the 
normalized  spatial  correlations  of  the  pressure  field  of  seven  different 
insonification  area  separations.  An  example  of  these  spatial  correlations 
for  surface  3  is  presented  in  Appendix  C.  The  surface  correlation  functions 
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300  0  300 

HEIGHTS  (cm) 


—  EPM  ~~  GAUSSIAN 


FIGURE  16 

THE  PROBABILITY  DENSITY  FUNCTION  FOR  THE  SECOND  ROUGH 
SURFACE  VIA  THE  EXTENDED  PRONY  METHOD 

STANDARD  DEVIATION  OF  THE  GAUSSIAN  CURVE  IS  215  cm 
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of  Figs.  18,  19,  and  20  are  monotonically  decreasing,  resembling 
exponential  functions.  Thus  an  exponential  least  squares  fit  between  the 
data  and  the  analytical  function  is  used  to  obtain  the  best  fitting  function. 
The  correlation  lengths  calculated  represent  the  e-folding  value  of  the 
exponential  curve.  The  correlation  function  fit  to  surface  1  is  reasonably 
close  and  yields  a  correlation  length  of  7  cm.  The  least  mean  square  fit  for 
surface  2  does  not  follow  the  structure  very  well,  and  if  an  exponential  is 
fit  visually  a  correlation  length  of  5.5  cm  is  obtained.  The  correlation 
length  for  surface  3  of  5.2  cm  is  also  reasonable. 

The  correlation  length  estimates  seem  to  agree  with  the  physical 
measurements  very  well.  Many  factors  have  contributed  to  this  agreement. 
The  pressure  field  is  neither  oversampled  nor  undersampled,  but  sampled 
sufficiently  to  yield  results  unadulterated  by  the  finite  extent  of  the  data. 
Also,  the  spatial  correlations  are  computed  for  spots  which  are 
incrementally  separated  a  fraction  of  a  correlation  length.  The  variation 
which  does  exist  is  believed  to  be  actually  due  to  the  variation  of  the 
correlation  along  different  orientations  (as  in  Fig.  3)  on  the  surface;  also, 
the  sampling  of  the  pressure  field  in  one  dimension  allows  calculation  of 
the  correlation  along  one  axis  of  the  surface. 

The  results  are  summarized  in  Table  I.  The  correlation  estimates  for 
surfaces  1,  2,  and  3  follow  very  closely  with  the  actual  physical 
measurements  of  the  correlation  lengths. 
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CORRELATION 
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SPOT  SEPARATION  (cm) 

*  DATA  ~~  EXPONENTIAL 


FIGURE  20 

CORRELATION  FUNCTION  FOR  THE  THIRD  ROUGH  SURFACE 
VIA  THE  POWER  SPECTRUM  OF  THE  SCATTERED  PRESSURE  FIELD 

THE  E-FOLDING  VALUE  OF  THE  EXPONENTIAL  CURVE  IS  5.18  cm 


VI.  CONCLUSIONS 


Several  techniques  have  been  outlined  which  have  yielded  predictions 
of  the  probability  density  function  of  surface  heights  and  correlation 
function.  The  estimates  of  these  statistical  characteristics  can  be 
obtained  directly  from  measurements  of  the  scattered  pressure  field  from 
the  rough  surfaces.  Thus  the  inverse  problem  has  been  investigated. 

The  correlation  function  estimates  agree  very  well  with  the 
measured  correlation  functions.  This  is  due  to  two  factors:  (I)  the 
insonification  spots  are  large  compared  to  the  dimensions  of  the 
correlation  length  and  (2)  the  pressure  field  spot  separation  and  spatial 
sampling  increment  are  both  a  fraction  of  a  correlation  length  in 
dimensioa  Therefore,  a  portion  of  the  surface  representative  of  the 
surface  correlation  function  has  been  insonified;  and  the  correlation 
function  estimate  has  been  observed  within  the  resolution  bounds 
necessary  to  determine  correlation  length.  Thus  Clay  and  lledwin's 
correlation  theory  was  verified,  this  time  with  stationary  randomly  rough 
surfaces. 

The  probability  density  estimates,  although  bounded,  were  not  as 
accurate  as  the  correlation  estimates.  The  estimated  rms  height  values 
were  all  exaggerated  with  the  exception  of  one  estimate.  From  the 
analysis  it  is  obvious  that  the  experimental  constraints  have  limited  the 


for  H,  the  complex  conjugate  transpose. 

Combining  Eqs.  B.I3  and  B.17  results  in 


q  q 

eji =  I  Z  (cgj-cgl!)rgHTh 

9=1  h=l 


such  that  0  will  have  one  eigenvector  Eq+1  which  is  orthogonal  to  the  q 
mode  vectors,  or 


B.20 


for  an  eigenvalue  of  zero. 

The  procedure  for  determination  of  q  is  to  fill  the  matrix  ©  to 
dimension  Q  by  Q  and  to  calculate  the  eigenvectors  and  eigenvalues.  If  L 
eigenvalues  are  equal  to  zero(or  equal  to  mo2  in  the  case  of  data  with  zero 
mean  noise  and  variance  a2),  then  q=Q-L.  The  matrix  could  be  recomputed 
to  order  q*1  by  q*1  and  the  eigenvectors  and  eigenvalues  regenerated.  From 
Eqs.  B.I6  and  B.20  it  is  noted  that  the  eigenvector  corresponding  to  a  zero 
eigenvalue  is  the  vector  with  the  coefficients  of  the  difference  equation  of 
Eq.  B.15.  Therefore  not  only  is  q  achieved  but  also  the  (aj  >. 

The  polynomial  equation  of  B.7  is  rooted  to  find  the  {zj }  and  Eq.  B.2 
reduces  to  a  set  of  linear  equations  which  can  be  written  in  matrix  form, 
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vector  Yh  starting  at  the  jth  £  step.  Recalling  Eq.  B A,  the  mode  vectors 
can  be  represented  with 


m  iT 


rh  =  li  V- V*] 


B.14 


for  which  there  are  q  vectors. 

Now,  seeking  a  solution  for  the  difference  equation  from  Eq.  B.IO, 


q 

1  Q(JL-j)  =  0 . 4=q,  ...  fl-1  .  B.15 

j=0 


A  pseudo- inverse  solution  is  achieved  for  which  Eq.  B.I5  takes  on  the  form 


0  A=  0 


B.I6 


where  ©  is  a  q+l  by  q+1  Hermitian  matrix  such  that 


ei»  =0iH°j 


B.I7 


and 


A  =  l  a0  a, ... 


B.I8 


rooted  to  obtain  the  {Zj}.  Subsequently,  the  {cj}  can  be  determined  by 

solving  the  equations  of  B.2  exactly  or  by  least  squares.  But  a  successful 
application  of  the  EP  method  depends  upon  knowledge  of  the  value  for  q. 
Many  methods  have  been  outlined  for  the  determination  of  q  and  the 
eigenanalysis  method  outlined  by  Van  Blaricum  and  Mittra  40  is  relatively 
easy  to  implement. 

We  begin  with  the  Dj  data  vectors 

Dj  =  [  Q(q-j)  Q(q-j+1)  ...  Q(q-j+m)  )T  B.11 

for  m+n-q-l  and  T  the  transpose.  These  can  also  be  written 

d 

Dj  =  Z  ^  explsh(q-j)A£]  exp[sh4A£)  ,  ft=0 . m  B.I2 

h=l 

or,  more  simply, 

d 

Dj  =  l  ChjTn  ,  B.B 

h=l 

where  Cbj=Cbexp[Sb(d_j)A^]  represents  the  coefficient  for  the  hth  mode 


such  that  a0  =  -1  and 


(z-z,)(z-z2)...  (z-z^)  =  0  .  B.8 

A  procedure  for  calculation  of  the  {ap  is  now  outlined.  The  first 

equation  in  B.2  (ft=0)  is  multiplied  by  a^,  the  second  by  aq_| .  and  the 

(q-H)th  by  a0=-1  and  the  sum  of  the  equations  is  computed.  Since  each  Zj 
satisfies  Eq.  B.8  the  result  is 


Q(q)  -  atQ(q-l)  -  ...  -  3^0(0)  =  0  B.9 

A  set  of  n-q  linear  equations  are  thus  obtained  by  using  thb  procedure  on 
the  remaining  equations  of  B.2.  The  resulting  equation  takes  on  the  form  of 
a  familiar  set  of  linear  prediction  equations, 

q 

QUO  =  2  3j  QU-j)  .  1  =  q . M-l  ,  B.10 

j=l 

which  can  be  solved  exactly  for  the  {aj },  if  M=2q,  or  approximately  by  a 
least  squares  estimation 

After  determination  of  the  {ap,  the  polynomial  equation  of  B.7  is 
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where  the  «j  is  a  damping  factor.  If  we  also  assume  that  the  {Cj}  are  also 
complex  amplitudes, 

Cj  =  Aj  exptiej]  ,  B.5 

then  the  FT  of  Eq.  B.2  is 

q 

<o(C)  =  1  Cj  2<xj/l«j2  4  {27T(C-Cj  )>2)  •  B.6 

j=» 

If  the  {Zj}  are  known,  then  Eq.  B.2  represents  a  set  of  M  linear 
equations  in  q  unknowns  to  be  solved  for  the  Cj.  For  M=q,  Eq.  B.2  can  be 

solved  exactly  as  Prony39  had  originally  intended,  and  for  M<q  a  linear 
least  squares  estimation  would  obtain  the  solution,  i.e.,  the  EP  method. 
Otherwise  the  determination  of  {Zj >  with  known  {Cj}  leads  to  the  solution 
of  a  set  of  nonlinear  equations. 

Consider  the  [z.)  to  be  the  roots  of  the  polynomial  equation 


APPENDIX  B 

DERIVATION  OF  THE  EXTENDED  PRONY  METHOD 


The  Prong  method  consists  of  expanding  the  data  set  known  at  evenly 
sampled  subintervals  with  a  basis  set  of  complex  exponentials,  or 

q 

Q(0  =  2  Cj  explSj$]  B.l 

j=1 

or.  for  discrete  values  of  £=AA£, 

q 

QUO  =  £  cj  [zj]1  .  A  =  0 . M-l  ,  B.2 

j=» 

( 

which  looks  similar  to  the  DFT  expansion.  However,  unlike  the  DFT  the  {sj} 
are  complex  in  general  and  non-harmonica! ly  related, 

Sj  =  «j  ♦  i27TCj  B.3 

and 

Zj=exptSjAO  ,  B.4 
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*(2j*1)  S(,,0j(2TT^C 

max-  ^max^ 


A.I8 


such  that  the  jth  Fourier  coefficient  of  the  solution  is 

<o>(0 .  j  >  =  <v(0 .  j  >/Xj  .  A.12 

Therefore,  a  solution  for  oi(0  is 

00 

o>(0=  1  t<v(0.'rtj>/Xj],nj(C)  A.13 

i=o 

for  "Cmax<C<+Cmax  and  the  series  is  convergent  in  the  mean  square  sense. 

The  j  are  related  to  the  prolate  angular  wave  functions  of  the  first 
kind,  Sm0j,  through 

^(0  =  s(t)oj^‘n’^mCmax.  C/CmaxV)ij  ,  A.14 

where 


=  ^  ^max  W2 


A.16 


<a ,  b>  =  a(0  b(0  dC  .  A.8 

J  -Cmax 

Although  not  explicitly  indicated,  the  <|>j  and  Xj  also  depend  upon  the 
product  2TT«mtmax. 

Slepian  and  Pollack38  describe  a  set  of  functions 

{T1|(0)  =  (/fxjT^(O)  A.9 

which  are  orthonormal  in  the  space  of  square  integrable  functions  over  the 
CL  range.  Multiplying  both  sides  of  Eq.  A.4  by  and  integrating  over 

“^max  to  +^max  9'ves 

♦Cmax 

<v(C).  t\ j(C)>  =  <  <p(x)sinc(C-x)dx.itj(C)>  A.io 

■  -Cmax 

The  symmetric  nature  of  the  kernel  and  Eq.  A.9  allows  Eq.  A.IO  to  be 
rewritten  as 


<v(C)  .  Tl;(C)>  =  X|<<1)(C)  ,  *n  j(C)> 


A.II 


v(C)  =  u>(x)  sinc(C-K)  dx  A.4 

J  -Cmax 

The  kernel  of  this  integral  equation  therefore  has  a  countably  infinite  set 
of  eigenfunctions 

{+j(C)).j  =  0,1 .  A.5 

corresponding  to  positive  eigenvalues 

X0  >  X,  >  X2  > ...  >  0  A.6 

such  that  the  set  {<J>j(C)J  is  complete  in  the  space  of  square  integrable 
functions  over  the  CL  range.  Also  the  are  orthogonal  over  the  CL  range, 

C+j .  $4  >  =  °.  j  *  A 

=  I.  j  =  *  A.7 


for  j,  4=0,1 .  and  where  the  inner  product  is  defined  with 


APPENDIX  A 

DERIVATION  OF  THE  BOJARSKY-LEWIS  METHOD 

The  Bojarsky-Lewis  method  is  developed  for  £L  functions.  Since  the 
data  are  finite  in  extent,  it  is  assumed  that  a  rectangular  window  U(0  is 
chosen  such  that 

U(0  =  I  .  £mjn  <  £  s*£max 

=  0 ,  otherwise  A.I 

and,  when  applied  to  the  characteristic  function,  yields 

V(0  =  U(0  Q(0  .  A.2 

Assuming  that  0  is  known  over  the  $L  range  and  that  V(0  has  an  inverse  FT, 
Eq.  A.2  when  inverse  transformed  and  the  convolution  theorum  invoked 
becomes 

♦  00  '• 

v(0  =  o)(k)  sinc(t-x)  dx  A.3 

.  -00 

for  v,  the  inverse  FT  of  V. 

For  a  real  surface  C  is  bounded  and  «>  in  Eq.  A.3  is  replaced  with 


reflection).  Also  Stanton34  has  successfully  implemented  a  pressure 
amplitude  analysis  in  an  inverse  type  problem.  Powell37  used  this  method 
in  an  analysis  of  the  author's  scattered  pressure  data. 

In  general,  the  objectives  of  the  inverse  problem  in  estimating  the 
statistical  parameters  of  a  randomly  rough  scattering  surface  were 
accomplished  with  some  limited  success. 
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resolutioa  Perhaps  the  most  critical  constraint  is  the  imposition  of  a 
very  small  beamwidth  source  to  insonify  small  spots  on  the  surfaces.  This 
results  in  a  scattered  pressure  field  which  is  very  limited  in  angular  extent 
resulting  in  a  very  small  observation  window.  With  the  additional 
constraint  of  a  single  frequency  of  incident  pressure,  this  limits  the 
observation  window  of  the  wave  number  space  of  the  characteristic 
function  of  surface  heights.  An  additional  limitation  to  the  resolution  is 
introduced  by  the  choice  of  near-normal  incidence  to  minimize  shadowing 
and  neglecting  surface  slopes  and  multiple  scattering. 

Various  improvements  can  be  made  to  increase  the  resolution  of  the 
experiment  and  as  a  result  the  accuracy  of  the  PDF  estimates.  The  first 
improvement  is  of  course  to  use  the  current  acoustic  model  and  simply 
increase  the  observation  window.  But  this  would  mean  using  a  larger 
beamwidth  source,  thus  requiring  a  larger  model  rough  surface  so  that  a 
representative  number  of  spots  could  be  insonified.  Secondly,  higher 
resolution  spectrum  analysis  techniques  could  be  used,  for  example, 
maximum  likelihood,  maximum  entropy,  etc. 

A  third  alternative  is  to  use  a  different  acoustic  scattering  model. 
This  could  be  done  by  modifying  the  present  theory  to  account  for 
broadband  incident  energy,  shadowing,  slopes,  or  multiple  scattering.  The 
result,  though,  would  be  most  certainly  non-trivial.  Clay  and  Medwin'9 
were  successful  in  applying  the  present  theory  for  finding  the  PDF  when 
they  used  various  frequencies  of  incident  pressure  to  sample  the  surface 
wave  number  space  (instead  of  varying  the  angles  of  incidence  and 
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APPENDIX  C 

SPATIAL  CORRELATION  FUNCTION 


The  spatial  correlation  function  of  the  pressure  field  contains 
information  which  will  allow  the  inference  of  the  correlation  function  of 
the  scattering  surface.  In  order  to  study  the  form  of  the  spatial 
correlation  function  suggested  by  Clay  and  Medwin,  the  spatial  pressure 
field  was  recorded  for  the  insonification  of  different  areas  on  the  model 
rough  surface.  The  pressure  field  was  then  spatially  correlated  by 
computing  the  correlation  between  the  pressure  fields  of  spots  which  have 
been  separated  a  given  distance  (for  all  possible  combinations  of  spots 
being  separated  that  distance).  Figure  C-1  represents  the  total  number  of 
spots  which  were  insonified  on  surface  3.  All  possible  combinations  of  the 
various  spot  separations  are  also  listed  by  the  number  of  separated  spot 
pairs.  Figures  C-2  through  C-8  represent  the  spatial  correlations  of  the 
pressure  fields  from  each  spot  separation.  It  should  be  noted  that  for 
greater  separations  the  maximum  of  the  correlation  function  decreases  and 
is  shifted  from  the  origin  an  amount  proportional  to  the  spot  separation. 


Spot 

Separation  (cm)  0.0  2.29  4.58  6.87  9.16  11.5  13.7 

*  Spot  Pairs  69  42  53  36  37  22  21 


figure  c-i 

LOCATION  OF  THE  69  SURFACE  INSONIFICATION  AREAS 
(-3  dB  SPOT  AREA)  ON  THE  THIRD  ROUGH  SURFACE 
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FIGURE  C-3 

SPATIAL  CORRELATION  OF  THE  SCATTERED  PRESSURE 
FIELD  FOR  SPOTS  SEPARATED  BY  2.29  cm 
THE  MAXIMUM  CORRELATION  VALUE  IS  0.695  AT  A  SPATIAL  LAG  OF  3.64  cm 
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SPATIAL 

CORRELATION 


-30.30  0.00  30.30 

SPATIAL  LAG  (cm) 


FIGURE  C-4 

SPATIAL  CORRELATION  FUNCTION  OF  THE  SCATTERED  PRESSURE 
FIELD  FOR  SPOTS  SEPARATED  BY  4.58  cm 

THE  MAXIMUM  CORRELATION  VALUE  IS  0.604  AT  A  SPATIAL  LAG  OF  6.06  cm 
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